#include "MerOverlapper.hh"


// Global buffer for performance in I/O. 
#define SEQ_BUFFER_SIZE 500000
char gseqbuf[SEQ_BUFFER_SIZE];


int bTable[MAX_BAND_SIZE][MAX_SEQ_LENGTH];

// Score for things that fall outside of the band.  
#define OUTSIDE_BAND_SCORE -100000

int gBandCount = 0;
//int mGapScore = -2;
//int mMatchScore = 1;
//int mMisMatchScore = -1;


//-------------------------------------------
// * ReadRepeatMers
//-------------------------------------------
//
//
void 
MerOverlapper::ReadRepeatMers(istream &fin,int minCount)
{
  long merCount = 1,addCount = 0;
  do{
    fin.getline(readbuf,READ_BUFFER_SIZE);
    // Find the first whitespace
    int i;
    for(i = 0;i< READ_BUFFER_SIZE;i++){		
      if ((readbuf[i] == '\t') | (readbuf[i] == ' ')) break;
    }

    // Convert remainder of the line to an int. 
    int count = atoi(readbuf+i+1);

    if (count >= minCount){
      mRepeatMers.Add(readbuf,count);
      addCount++;
    }
    if (!(merCount%500000)){
      cout<<"Mers seen:"<<merCount;
      cout<<"\tInsertions: "<<mTotalInsertions;
      cout<<"\tCollisions: "<<mCollisions;
      double CollisionRate=(double)mCollisions/(double)mTotalInsertions;
      cout<<" ("<<CollisionRate<<")";
		
      cout<<"\tWalk steps: "<<mWalkSteps;
      double AvgWalkSteps = (double)mWalkSteps/(double)mCollisions;
      cout<<" ("<<AvgWalkSteps<<")";

      double loadFactor = (double)mTotalInsertions/(double)HASHSIZE;
      cout<<"\tLoad: "<<loadFactor<<endl;
    }
    merCount++;

  }while(!fin.eof());
}

//-------------------------------------------
// * ComputeAllOverlaps
//-------------------------------------------
// Compare each read in the file qin to the MerHash 
// analyzing the common mers between each read and 
// the reads in the hash. The overlap results are 
// written as a sparse graph to fout.  
// 
void 
MerOverlapper::ComputeAllOverlaps(istream &qin,ostream &fout)
{
  // Read each read from qin, and determine it's overlaps 
  // with sources in the hash, writing results to fout. 
  cout<<"\tComparing reads to source..."<<endl;
  int readCount=0;
  bool readGood;
  gseqbuf[0]='1';// Denotes first read of a file. 
  do{
    Sequence3Bit QueryRead;	 
    readGood = QueryRead.GetNextReadTrimLowQual(qin,gseqbuf,SEQ_BUFFER_SIZE);
    if ((readGood) && (QueryRead.mBaseSize > 0)){		
      FindOverlapsWithSource(QueryRead,fout);		
      readCount++;
      if (!(readCount%25000)){
	cout<<"\t"<<readCount<<" reads examined."<<endl;
      }
      // cout << "26:  Read " << readCount << " passed!" << endl;
    }
    //else{
    //  cout << "26:  Read " << readCount << " failed!" << endl;
    //}
  }while(!qin.eof() && readGood);
  cout<<"\tdone. ReadCount="<<readCount<<endl;
  cout<<"\t"<<gBandCount<<" banded alignments performed."<<endl;
}

//-------------------------------------------
// * SaveOverlaps
//-------------------------------------------
//
void
MerOverlapper::SaveOverlaps(ostream &fout,string &QueryName,int ReadSize,ReadOverlaps &Overlaps)
{
  if (Overlaps.size() == 0) return;
  bool firstRead = true;

  ReadOverlaps::iterator i;
  vector<string> &names = mSourceHash.mSourceNames;
  string sourceName;
  for(i = Overlaps.begin();i != Overlaps.end();i++){
    int seqIdx = (*i).first;
    OverlapInfo eInfo = (*i).second;
    sourceName = names[seqIdx];
	 
    if (mFlipDirection) {
      fout << sourceName << "\t" << QueryName << " {" <<
	eInfo.mSpan << "," << eInfo.mScore << ",";
      if (mOriginRelative) {
	fout <<eInfo.mLeftExtension<<","<<eInfo.mRightExtension;
      }
      else {
	if (eInfo.mSense == 'f') {
	  fout << -eInfo.mLeftExtension << "," << -eInfo.mRightExtension;
	}
	else {
	  fout << -eInfo.mRightExtension << "," << -eInfo.mLeftExtension;
	}
      }
      fout<<","<<eInfo.mSense<<","<<eInfo.mCopyNum<<"}" << endl;
    }
    else {
      // Write source name if first read. 
      if (firstRead){
	fout<<QueryName<<"\t";
	firstRead = false;
      }

      fout<<sourceName<<" {"<<eInfo.mSpan<<","<<eInfo.mScore<<",";
      if (mOriginRelative) {
	if (eInfo.mSense == 'f') {
	  fout << -eInfo.mLeftExtension << "," << -eInfo.mRightExtension;
	}
	else {
	  fout << -eInfo.mRightExtension << "," << -eInfo.mLeftExtension;
	}
      }
      else {
	fout <<eInfo.mLeftExtension<<","<<eInfo.mRightExtension;
      }
      fout<<","<<eInfo.mSense<<","<<eInfo.mCopyNum<<"}\t";
    }
  }

  // Don't want to write a newline unless something was actually saved. 
  if (!firstRead) fout<<endl;   
}


//-------------------------------------------
// * FindOverlapsWithSource
//-------------------------------------------
// Finds the overlaps between a single query and 
// the reads in the hash. 
// 
void
MerOverlapper::FindOverlapsWithSource(Sequence3Bit &QueryRead,ostream &fout)
{
  // Skip trivial reads. 
  if (QueryRead.mBaseSize < mMinReadSize){
    cout<<QueryRead.mName<<" size < "<<mMinReadSize<<" skipping."<<endl;
    return;
  }

  MerSearchResultVec<MerLocation> FwdHits,RevHits;

  // Choose sampling method
  if (mInteriorSampling){
    //cout << "26:  Searching, method 1..." << endl;
    SearchSamplesAgainstHash(QueryRead,FwdHits,RevHits); 
  }else{
    // Sample QueryRead, and search samples against hash.
    //cout << "26:  Searching, method 2..." << endl;
    SearchEndSamplesAgainstHash(QueryRead,FwdHits,RevHits); 
  }

  // Sort hits.  KJD This could be made faster.
 
  // Each sort is stable, and preserves the sub-order of the 
  // previous sort.  So, source read being equal, it will 
  // sort by the MerCount, and MerCount being equal it will 
  // reverse sort by the base idx (i.e. want to favor an 
  // interior hit as a seed over an end hit.  
  FwdHits.SortBySourceBaseIdx();
  FwdHits.SortByMerCount();
  FwdHits.SortBySourceRead();

  RevHits.SortBySourceBaseIdx();
  RevHits.SortByMerCount();
  RevHits.SortBySourceRead();



  //  cout<<endl;
  //cout<<"==========================="<<endl;
  //FwdHits.Print(cout);
  //cout<<"====== Reverse hits ======"<<endl;
  //RevHits.Print(cout);
  //cout<<"==========================="<<endl;


  // Analyze hits to determine overlaps.  
  ReadOverlaps Overlaps;

  EvaluateFwdHits(FwdHits,Overlaps,QueryRead);
  EvaluateRevHits(RevHits,Overlaps,QueryRead);

  // Write overlaps to disk. 
  SaveOverlaps(fout,QueryRead.mName,QueryRead.mBaseSize,Overlaps);
}


void 
MerOverlapper::InitGapPenalties(int band)
{	 	 
  int btrow,btcol;
  register int NumBandRows = 2*band+1;

  // first clear everything
  for(btrow = 0;btrow < MAX_BAND_SIZE;btrow++){
    for(btcol = 0;btcol<  MAX_SEQ_LENGTH;btcol++){
      bTable[btrow][btcol] = 0;
    }
  }

  // Init table with gap penalties. 
  for(btrow = band;btrow < NumBandRows;btrow++){
    bTable[btrow][0] = mGapScore*(btrow - band);
  }
  
  for(btcol =0;btcol <= band;btcol++){
    bTable[band-btcol][btcol] = btcol * mGapScore;
  }
}

void
MerOverlapper::Bug2()
{
  cout<<"You've found a bug in BandedAlignment"<<endl;
}


//-------------------------------------------
// * ComputeBandedAlignmentScore
//-------------------------------------------
// This does a banded alignment without keeping a traceback. 
// This really should use the version in Aligner.hh, but 
// before I can do that I have to unify the Sequence3Bit and 
// the Sequence interface (or, pending performance tests, 
// replace Sequence entirely with Sequence3Bit).  So, for now, 
// this rather awkwardly lives here....
// 
void 
MerOverlapper::ComputeBandedAlignmentScore(Sequence3Bit &Query,int QStart,int QEnd, 
					   Sequence3Bit &Source,int SStart,int SEnd, int band, int &score)
{

  gBandCount++;

  register int btcol,btrow;
  register int diagScore,leftScore,topScore;
  register int NumBandRows = 2*band+1;
  register int lastBandRow =2*band;
  register int seqSize = QEnd-QStart+1;
  
  // A banded alignment looks at only a subset of the alignment table 
  // for performance.  In the code that follows, the prefix bt stands for 
  // band table.  All variables with bt prefix are in banded alignment table 
  // coordinates.   There is an imaginary full alignment table, with it's own 
  // coordinates.  I call this the virtual table, and variables in that coordinate
  // system are prefixed with "virtual". 

  register int rowStart = band;
  for(btcol = 1; btcol <= seqSize;btcol++){
    register int btcolm = btcol -1;

    for(btrow = rowStart; btrow < NumBandRows;btrow++){
      register int virtualRow = btcol - band + btrow;

      register int btrowm = btrow -1;
      register int btrowp = btrow +1;

      // Offset into sequence. 
      register int basecol = btcol+QStart-1;
      register int baserow = virtualRow+SStart-1;

      // There is no score in the lower band size triangle of the matrix
      // Trying to access the Source there will, if the SEnd is the last character in the 
      // source, go past the end of the source.  An equally valid test would be to test
      // btrow and btcol imply that you are in the lower right band sized triangle of the 
      // band table. 
      if (baserow < Source.mBaseSize){
	if (Query[basecol] == Source[baserow]){
	  diagScore = bTable[btrow][btcolm]+mMatchScore;
	}else{
	  diagScore = bTable[btrow][btcolm]+mMisMatchScore;
	}
		  
	if (btrow != lastBandRow) leftScore = bTable[btrowp][btcolm]+mGapScore;
	else leftScore = OUTSIDE_BAND_SCORE; 
		  
	if (virtualRow > 0) topScore = bTable[btrowm][btcol]+mGapScore;
	else topScore = OUTSIDE_BAND_SCORE;

		
	// Find the max of the scores and record.  
	if (diagScore > leftScore){
	  if (diagScore > topScore){
	    bTable[btrow][btcol]=diagScore;
	  }else{
	    bTable[btrow][btcol]=topScore;
	  }
	}else{
	  if (leftScore > topScore){
	    bTable[btrow][btcol] = leftScore;
	  }else{
	    bTable[btrow][btcol] = topScore;
	  }
	}
      }
    }
	 
    rowStart--;
    if (rowStart < 1) rowStart = 1;

  }
  score = bTable[band][seqSize];  
}

//-------------------------------------------
// * EvaluateFwdHits
//-------------------------------------------
// Uses the first hit for each source-query pair as an anchor to compute the boundaries
// for a banded alignment.  Performs the banded alignment, and computes the left and right
// extensions.    A negative left or right extension means that the source can not be 
// extended in that direction, but is reported anyway so that the magnitude of the overhang
// in the other direction can be known.  A positive left or right extension means that the 
// source can be extend in the left or right direction by that ammount.  
//
//
void
MerOverlapper::EvaluateFwdHits(MerSearchResultVec<MerLocation> &Results,ReadOverlaps &Overlaps,
			       Sequence3Bit &Query)
{
  if (Results.size() == 0) return;

  int SampleSpacing = mSampleGap;
  //int closeToBeginning = mSamplesPerEnd * SampleSpacing;  

  int r = 0;
  int rsize = Results.size();
  do{

    // For a given source Idx, want to find the hit that is based on 
    // the rarest n-mer, out of hits for a particular source/query pair. 

    register int CurrentSourceIdx = Results[r].SourceRead();	 
    int SourceStart = Results[r].SourceBaseIdx();
    int QueryStart = Results[r].mQueryBaseIdx;

    int rcount = Results[r].mMerGlobalFrequency;

    do{
      r++;
      if (r >= rsize) break; // End this loop if we have exhausted results
    }while(Results[r].SourceRead() == CurrentSourceIdx);	
	 
    Sequence3Bit &Source = mSourceHash.mSourceSeqs[CurrentSourceIdx];
    int score=-10000;


    // The amount of the source and query that extends to the right of the sample
    int SourceRight, QueryRight;
    // The degree to which the source can be extended to the left or right by this overlap. 
    int LeftExtension, RightExtension;
    // Coordinates for the banded alignment. 
    int QAlignStart,QAlignEnd,SAlignStart,SAlignEnd;
	 
    SourceRight = Source.mBaseSize - SourceStart;
    QueryRight = Query.mBaseSize - QueryStart;

    // There are four cases to attend to.  
    if ((QueryRight>SourceRight) && (SourceStart>QueryStart)){
      // Case I
      // S:  ---------->
      // Q:       ---------->
      QAlignStart = 0;
      SAlignStart = SourceStart - QueryStart;
      SAlignEnd = Source.mBaseSize - 1;
      QAlignEnd = SAlignEnd - SAlignStart;

      LeftExtension = -SAlignStart;
      RightExtension = Query.mBaseSize - QAlignEnd -1;

    }else if ((SourceRight>QueryRight) && (QueryStart > SourceStart)){
      // Case II
      // S:       ---------->
      // Q:  ---------->
      SAlignStart = 0;
      QAlignStart = QueryStart - SourceStart;
      QAlignEnd = Query.mBaseSize-1;
      SAlignEnd = QAlignEnd - QAlignStart;

      LeftExtension = QAlignStart;
      RightExtension = -(Source.mBaseSize - SAlignEnd -1);

    }else if ((SourceRight>=QueryRight) && (SourceStart>=QueryStart)){
      // Case III
      // S:  ---------->
      // Q:   ----->
      QAlignStart = 0;
      SAlignStart = SourceStart - QueryStart;
      QAlignEnd = Query.mBaseSize - 1;
      SAlignEnd = SourceStart + (QAlignEnd - QueryStart);

      LeftExtension = -SAlignStart;
      RightExtension = -(Source.mBaseSize - SAlignEnd -1);

    }else if ((QueryRight>=SourceRight) && (QueryStart>=SourceStart)){
      // Case IV
      // S:    ----->
      // Q:  ---------->
      SAlignStart = 0;
      QAlignStart = QueryStart - SourceStart;
      SAlignEnd = Source.mBaseSize - 1;
      QAlignEnd = QueryStart + (SAlignEnd - SourceStart);

      LeftExtension = QAlignStart;
      RightExtension = (Query.mBaseSize - QAlignEnd -1);

    }else{
      cout<<"Missing case:  QueryRight="<<QueryRight<<" SourceRight="<<SourceRight;
      cout<<" SourceStart="<<SourceStart<<" QueryStart="<<QueryStart<<endl;
    }

    ComputeBandedAlignmentScore(Query,QAlignStart,QAlignEnd,Source,SAlignStart,SAlignEnd,mBandSize,score);

    int span = QAlignEnd - QAlignStart+1;	 
    int cutoff = (int) ((1 - ((float)2.0*mMaxMismatchPercent/100.0)) * span);	
    if (score >= cutoff){

      // If no such edge exists previously, this must be the best one found so far so add it
      // This is so that both fwd and reverse samples have a crack at being the best.   
      int escore = Overlaps[CurrentSourceIdx].mScore;
      if (escore == 0){		  
	Overlaps[CurrentSourceIdx].mScore = score;
	Overlaps[CurrentSourceIdx].mSpan = span;
	Overlaps[CurrentSourceIdx].mLeftExtension = LeftExtension;
	Overlaps[CurrentSourceIdx].mRightExtension = RightExtension;
	Overlaps[CurrentSourceIdx].mSense = 'f';
	Overlaps[CurrentSourceIdx].mCopyNum = rcount;
      }else{
	// If the new score is better
	if (escore < score){
	  Overlaps[CurrentSourceIdx].mScore = score;
	  Overlaps[CurrentSourceIdx].mSpan = span;
	  Overlaps[CurrentSourceIdx].mLeftExtension = LeftExtension;
	  Overlaps[CurrentSourceIdx].mRightExtension = RightExtension;
	  Overlaps[CurrentSourceIdx].mSense = 'f';
	  Overlaps[CurrentSourceIdx].mCopyNum = rcount;
	}
      }
    }	 
  }while(r < rsize);
}


//-------------------------------------------
// * EvaluateRevHits
//-------------------------------------------
//
void
MerOverlapper::EvaluateRevHits(MerSearchResultVec<MerLocation> &Results,ReadOverlaps &Overlaps,
			       Sequence3Bit &Query)
{
  if (Results.size() == 0) return;

  int SampleSpacing = mSampleGap;
  //int closeToBeginning = mSamplesPerEnd * SampleSpacing;  

  int r = 0;
  int rsize = Results.size();
  do{

    register int CurrentSourceIdx = Results[r].SourceRead();	 
    int SourceStart = Results[r].SourceBaseIdx();
    int QueryStart = Results[r].mQueryBaseIdx;
    int rcount = Results[r].mMerGlobalFrequency;

    do{
      r++;
      if (r >= rsize) break; // End this loop if we have exhausted results
    }while(Results[r].SourceRead() == CurrentSourceIdx);	
	 
    Sequence3Bit &Source = mSourceHash.mSourceSeqs[CurrentSourceIdx];
    int score = -10000;
	 
    Sequence3Bit QueryRev;
    Query.ReverseComplement(QueryRev);

    // Query is Reverse so need to transform coordinates
    QueryStart = Query.mBaseSize - QueryStart - Mer::mMerSize;
	 	 
    // The amount of the source and query that extends to the right of the sample
    int SourceRight, QueryRight;
    // The degree to which the source can be extended to the left or right by this overlap. 
    int LeftExtension, RightExtension;
    // Coordinates for the banded alignment. 
    int QAlignStart,QAlignEnd,SAlignStart,SAlignEnd;

    SourceRight = Source.mBaseSize - SourceStart;
    QueryRight = Query.mBaseSize - QueryStart;

    // There are four cases to attend to.  
    if ((QueryRight>SourceRight) && (SourceStart>QueryStart)){
      // Case I
      // S:  ---------->
      // Q:       <----------
      QAlignStart = 0;
      SAlignStart = SourceStart - QueryStart;
      SAlignEnd = Source.mBaseSize - 1;
      QAlignEnd = SAlignEnd - SAlignStart;

      LeftExtension = -SAlignStart;
      RightExtension = Query.mBaseSize - QAlignEnd -1;

    }else if ((SourceRight>QueryRight) && (QueryStart > SourceStart)){
      // Case II
      // S:       ---------->
      // Q:  <----------
      SAlignStart = 0;
      QAlignStart = QueryStart - SourceStart;
      QAlignEnd = Query.mBaseSize-1;
      SAlignEnd = QAlignEnd - QAlignStart;

      LeftExtension = QAlignStart;
      RightExtension = -(Source.mBaseSize - SAlignEnd -1);

    }else if ((SourceRight>=QueryRight) && (SourceStart>=QueryStart)){
      // Case III
      // S:  ---------->
      // Q:   <-----
      QAlignStart = 0;
      SAlignStart = SourceStart - QueryStart;
      QAlignEnd = Query.mBaseSize - 1;
      SAlignEnd = SourceStart + (QAlignEnd - QueryStart);

      LeftExtension = -SAlignStart;
      RightExtension = -(Source.mBaseSize - SAlignEnd -1);

    }else if ((QueryRight>=SourceRight) && (QueryStart>=SourceStart)){
      // Case IV
      // S:    ----->
      // Q:  <----------
      SAlignStart = 0;
      QAlignStart = QueryStart - SourceStart;
      SAlignEnd = Source.mBaseSize - 1;
      QAlignEnd = QueryStart + (SAlignEnd - SourceStart);

      LeftExtension = QAlignStart;
      RightExtension = (Query.mBaseSize - QAlignEnd -1);

    }else{
      cout<<"Missing case:  QueryRight="<<QueryRight<<" SourceRight="<<SourceRight;
      cout<<" SourceStart="<<SourceStart<<" QueryStart="<<QueryStart<<endl;
    }

    ComputeBandedAlignmentScore(QueryRev,QAlignStart,QAlignEnd,Source,SAlignStart,SAlignEnd,mBandSize,score);

    int span = QAlignEnd - QAlignStart+1;	
    int cutoff =  (int) ((1 - ((float)2.0*mMaxMismatchPercent/100.0)) * span);	
    if (score >= cutoff){

      // If no such edge exists previously, this must be the best one found so far so add it
      // This is so that both fwd and reverse samples have a crack at being the best. 
      int escore = Overlaps[CurrentSourceIdx].mScore;
      if (escore == 0){		  
	Overlaps[CurrentSourceIdx].mScore = score;
	Overlaps[CurrentSourceIdx].mSpan = span;
	Overlaps[CurrentSourceIdx].mLeftExtension = LeftExtension;
	Overlaps[CurrentSourceIdx].mRightExtension = RightExtension;
	Overlaps[CurrentSourceIdx].mSense = 'r';
	Overlaps[CurrentSourceIdx].mCopyNum = rcount;
      }else{
	// If the new score is better
	if (escore < score){
	  Overlaps[CurrentSourceIdx].mScore = score;
	  Overlaps[CurrentSourceIdx].mSpan = span;
	  Overlaps[CurrentSourceIdx].mLeftExtension = LeftExtension;
	  Overlaps[CurrentSourceIdx].mRightExtension = RightExtension;
	  Overlaps[CurrentSourceIdx].mSense = 'r';
	  Overlaps[CurrentSourceIdx].mCopyNum = rcount;
	}
      }  
    }
  }while(r < rsize);
}


//-------------------------------------------
// * SearchSamplesAgainstHash
//-------------------------------------------
// Sample QueryRead, and search samples against hash.  
// Searches both forward and reverse samples.  Skips samples
// that contain N or X.  Samples are taken from left and right 
// ends. 
//
void
MerOverlapper::SearchSamplesAgainstHash(Sequence3Bit &QueryRead,
					MerSearchResultVec<MerLocation> &FwdLocations,
					MerSearchResultVec<MerLocation> &RevLocations)
{  
  // Search samples from the left end of read against hash.  
  // Find the first all-base (no N or X) mer
  int merStart = 0;
  int merSize = Mer::mMerSize;
  int querySize = QueryRead.mBaseSize;

  while ((merStart + merSize) <= querySize &&
	 QueryRead.ContainsNX(merStart, merSize)) 
    {
      merStart += mSampleGap;
    }
  
  //cout<<"Initial merStart="<<merStart<<endl;
  for(int i = 0;i< mSamplesPerRead;i++){
    if ((merStart + merSize) <= querySize){
      mSourceHash.FindMerLocations(QueryRead,merStart,FwdLocations,mRepeatMers);
      mSourceHash.FindRevMerLocations(QueryRead,merStart,RevLocations,mRepeatMers);
    }

    // Find next non NX sample
    do {
      merStart += mSampleGap;
    } while((merStart + merSize) <= querySize && 
	    QueryRead.ContainsNX(merStart, merSize));
  }
}


//-------------------------------------------
// * SearchEndSamplesAgainstHash
//-------------------------------------------
// Sample QueryRead, and search samples against hash.  
// Searches both forward and reverse samples.  Skips samples
// that contain N or X.  Samples are taken from left and right 
// ends. 
//
void
MerOverlapper::SearchEndSamplesAgainstHash(Sequence3Bit &QueryRead,
					   MerSearchResultVec<MerLocation> &FwdLocations,
					   MerSearchResultVec<MerLocation> &RevLocations)
{  
  int SamplesPerEnd = (int) (mSamplesPerRead/2.0);

  // Search samples from the left end of read against hash.  
  // Find the first all-base (no N or X) mer
  int merStart = 0;
  int querySize = QueryRead.mBaseSize;
  int merSize = Mer::mMerSize;

  while ((merStart + merSize) <= querySize &&
	 QueryRead.ContainsNX(merStart, merSize)) 
    {
      merStart += mSampleGap;
    }
  
  //cout <<"Initial merStart="<<merStart<<endl;
  //cout << "26:  Starting at " << merStart << "..." << endl;
  int leftMerEnd;
  for(int i = 0;i< SamplesPerEnd;i++){
    if ((merStart + merSize) <= querySize){
      //cout << "26:  Finding forward locations..." << endl;
      mSourceHash.FindMerLocations(QueryRead,merStart,FwdLocations,mRepeatMers);
      //cout << "26:  Finding reverse locations..." << endl;
      mSourceHash.FindRevMerLocations(QueryRead,merStart,RevLocations,mRepeatMers);
    }
    //else{
    //  cout << "26:  The query size " << querySize << " was too small." << endl;
    //}
    // Endpoint of last attempted mer from the left. 
    leftMerEnd = merStart + merSize - 1;
	 
    do{
      merStart += mSampleGap;
    } while((merStart + merSize) <= querySize && 
	    QueryRead.ContainsNX(merStart, merSize));
  }
	 
  // Search samples from the right end of read against hash.
  // Find the first all-base mer. 
  merStart = querySize - merSize;
  while ((merStart > leftMerEnd) && 
	 QueryRead.ContainsNX(merStart, merSize))
    {
      merStart -= mSampleGap;
    }
  //cout<<"Right merStart="<<merStart<<endl;

  for(int i = 0;i< SamplesPerEnd;i++){
    if (merStart > leftMerEnd){
      mSourceHash.FindMerLocations(QueryRead,merStart,FwdLocations,mRepeatMers);
      mSourceHash.FindRevMerLocations(QueryRead,merStart,RevLocations,mRepeatMers);
    }
		
    do{
      merStart -= mSampleGap;
    } while((merStart > leftMerEnd) && 
	    QueryRead.ContainsNX(merStart, merSize));
  }
}
